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Abstract. A stochastic simulation algorithm for the computation of multitimc 
correlation functions which is based on the quantum state diffusion model of 
open systems is developed. The crucial point of the proposed scheme is a suitable 
extension of the quantum master equation to a doubled Hilbert space which is 
then unraveled by a stochastic differential equation. 
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Within the framework of the recently developed stochastic wave function approach 
to open quantum systems g, [| [|, [|, ||, [j], ^] the state of a system is not described 
by a reduced density matrix but by a pure stochastic state vector tp t whose covariance 
matrix is equal to the reduced density matrix of the system. One of these models 
which was motivated by a dynamical description of the measurement process || is 
the quantum state diffusion model introduced by Gisin and Percival ^ [j). In this 
approach, the time-evolution of the wave function ipt is governed by the Ito stochastic 
differential equation 

#t = - \H4, t dt +\Y. [ 2 ( L ]Ut L 3 - L ) L J - fo>*t ^tdt 
j 

+ £fo-<-&i>*]lM&t. (!) 
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where {Lj)^ t is a short-hand notation for (ip t \L\ip t ), and d^jt is the differential of a 
complex valued Wiener process with means and correlations 

{dtit) = (d&d&t) = 0, (d^* t ) = Sijdt. (2) 

The operators H and Lj acting in the Hilbert space TC of the system are the free 
Hamiltonian and the Lindblad operators describing dissipation, respectively. The 
link to the density matrix description of open quantum systems is established - as 
mentioned above - through the covariance matrix of the stochastic wave function ipt, 
i. e., 

p t = E(\i> t ){ip t \). (3) 

The symbol E denotes the expectation value with respect of the stochastic processes 
ipt . The equation of motion of the density matrix p t is obtained by inserting eq. ([l]) 
in eq. (|^) which yields the quantum master equation 

P (t) = -i [H, p(t)] + ±J2 [ 2L 3p(t)L] ~ L]L jP (t) - P (t)L%] . (4) 
j 

This equation - or alternatively the stochastic differential equation (Q) - determines 
the time-evolution of one-time expectation values of system obscrvablcs. Multitimc 
correlation functions which are of special interest in quantum optics or in solid state 
physics are not specified by these equations. In order to define these quantities, wc 
will first define the matrix elements of some system operator A in the Heisenberg 
picture. In the density matrix approach, these matrix elements are defined through 
the quantum regression theorem jl0|, |ll| as 

A t (^ O) Vo) = (0o ) to|A(t)|Vo ) to)=Tr{^(i,i o ){|V'o)(0o|}}, (5) 

where V(t,to) is the time-evolution superoperator corresponding to the quantum 
master equation (^). Unfortunately the quantum regression theorem cannot be 
applied directly to the stochastic wave function approach, since the initial "density 
matrix" |?/;o)(</>o| is not necessarily Hermitian, and hence it can in general not be the 
covariance matrix of some stochastic wave function. This problem can be resolved by 
extending the quantum master equation into a doubled Hilbert space 7i — Ti. H. in 
the following way: we define a density matrix p(t) as 
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where pij (t) are operators on Ti. and accordingly replace the Hamiltonian H and the 
Lindblad operators Lj by the operators 



H = 



H 




Q 
H 



Lj 






u 



(7) 



in the doubled Hilbert space Ti.. Then we formulate the extended quantum master 
equation 



P^) = -i [H, Rtj\ [^jp{t)L\ - L}Ljp(t) - p{t)L)Lj 



(8) 



The crucial point of this construction is that each element £ij(t) of the density matrix 
p{t) is a solution of the original quantum master equation (|4]) . Consider now the initial 
condition 



P(t ) 



\Oo)(0o\ = ^ 



\<h)(<h\ \<h)(ih\ 



(9) 



where 6> = (<j>o, V'o) T /v / 2 is an element of the doubled Hilbert space Tt. (Throughout 
this letter the superscript T denotes the transpose). Obviously, the matrix elements 
of some operator A are then given by 



A t (0 o ,Vo) = 2Tr{Ap 21 (t)}. 



(10) 



By construction, the initial density matrix p(to) — \9o)(9q\ is positive and we may 
choose any unraveling of the extended quantum master equation (|^) by a stochastic 
process for the calculation of its time-evolution and hence for the calculation of 
operators in the Heisenberg picture (A similar idea has been proposed in Ref. Q 
Appendix D). 

Applying the above procedure to the quantum state diffusion model we obtain 
for example the equation of motion for the wave function 9 t — (<ftt,ipt) T £ TL in the 
Ito form 



d0 t = - iH9 t dt 



2 ^ 



j 



[Z - (L 3 ) e ] e t di 



e t dt 



(ii) 



The matrix elements of A are simply obtained as 

A t (<j> ,Tp ) = 2E eo ((<j> t \A\i; t )), (12) 
where E# denotes the expectation value with respect to the initial condition 6q. Note, 



that eq. (11) is constructed in such a way that the norm of the state vector 9 t is 
preserved, i. e., ||0t|| 2 = ||0t|| 2 + HV'tll 2 = 1- From a numerical point of view it is 
more efficient to drop this restriction and to work with unnormalized state vectors 9t , 
whose time-evolution is governed by the quasi-linear stochastic differential equation 



d§ t = -iH8 t dt + J2 L A (d£j t + (L]) 9t dt) - i ^ L)L 3 6 t dt. (13) 

3 3 

Accordingly, the matrix elements of the operator A are defined as 



A. 



(4,^o) = 2£ 9o ({$ t \A\$ t )/\\e t \f). 



(14) 
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As a particular example, we consider a two-level system with H = coupled to 
the vacuum using the Lindblad operator a~ , and calculate the matrix element 
(0o|c + (i)|7/>o), where 4> = (1,0) T and ip = (1, l) T /v / 2- The analytical solution 

(0ok + (Wo) = ^ e - t/2 (15) 

is readily obtained by integrating the quantum master equation. In Fig.|]we compare 
the numerical solution obtained using the scheme described above for 10 realizations 
(diamonds) with the analytical solution (thick line). Obviously, both solutions are in 
excellent agreement. 

Alternatively, Gisin proposes in Rcf. [jl3| a similar scheme for the calculation 
of matrix elements which is based on the coupled system of stochastic differential 
equations 



dip t = -iHiptdt + - [2lj(ip t ,4>t)*Lj - L^Lj - lj((t>t,il)t)lj{ipt,<t>t) 
j 

+ ^2[L - lj((f> t ,ipt)]iptd^jt, 
j 



ip t dt 



+ 5^[£-ii(lk,&)]&d&, (16) 

j 

where lj(a,P) = (a\Lj\/3)/(a\(3). These equations are constructed in such a way that 
the scalar product ((j>t\ipt) remains constant during the time-evolution of the system, 
i. e., the matrix element of the unity operator It — I are calculated correctly for each 
realization of the stochastic process (and not only in the mean). In addition, he also 
proposes in Ref. [|l3j a pair of quasi- linear equations, which could be used for the 
numerical simulation. However, although the above equations correctly reproduce the 
equation of motion for the matrix elements, the numerical integration of the stochastic 
differential equations for the system described above, suggests that these equations 
are not stable in general. In order to demonstrate this, we have also plotted in Fig. [t] 
the numerical solution of the quasi-linear stochastic differential equations for various 
step-sizes (h = 0.01, 0.001, 0.0001) and 10 4 realization each. The systematic deviation 
of the numerical and analytical solutions for t > 0.37" 1 is evident. We believe, that 
these deviations are due to the fact, that the solution of the deterministic part of 
the stochastic differential equation is unstable for this particular model which leads 
to immense fluctuations in the solution of the stochastic differential equation. Note, 
that the fluctuations are even much larger for the integration of the "unity- preserving" 
equation (|l^). 

The simulation algorithm in the doubled Hilbert space for the calculation of 
matrix elements in the Heisenberg picture is the basis for the computation of multitimc 
correlation functions such as g(t,t + r) = (t[>o\A(t + T)B(t)\ijjo) and we propose 
the following procedure: start in the state ipo and propagate it up to the time t 
using the stochastic differential equation (|l|) to obtain iji t . Define the state vector 
&t = C0t> Bip t ) T /yl + H-B^tH 2 and propagate it up to the time t + r by integrating 
the extended stochastic differential equation (pi]). The two-time correlation function 
g(t, t + t) is then given by 

9 (t, t + T) = E [(1 + HB^II 2 ) {'Pt+rWt+r)] ■ (17) 
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As a specific example we have computed the first order correlation function (er + (t + 
t)<t~ (t)) s for a coherently driven two-level atom on resonance in the steady state with 
Rabi frequency Q = IO7. To this end, we started with a random initial state vector ip 
drawn from a uniform distribution on 7i and propagated it up to t = 307 -1 in order 
to reach the steady state regime. Then we proceeded as described above. The result 
of the numerical simulation is shown in Fig. |^ (a) for 10 4 realizations. The numerical 
performance of the algorithm is demonstrated in Fig. ^ (b) where we have plotted 
the computational time which is necessary to obtain a given accuracy measured by 
the relative mean square error (solid line) and the estimated standard deviation of the 
samples (dashed line) . These results are compared with an alternative procedure which 
is based on an unraveling of the extended quantum master equation by a piecewise 
deterministic jump process (see Ref. |l2) ). The algorithm based on quantum jumps is 
about two times faster than the one based on the quantum state diffusion model. At 
a first glance, this result is surprising, since the individual realizations of the diffusion 
process are smooth and "closer" to the real solution. But this is outweighed by the 
fact that for the integration of the stochastic differential equation we have to draw 
two random numbers per time step and Lindblad operator, whereas in the quantum 
jump method we have to generate only two random numbers per jump. Thus, a single 
realization of the diffusion process is more accurate, but takes longer to be computed. 

To summarize, we have shown that operators in the Heisenberg picture and 
multitime correlation functions can be calculated within the framework of the quantum 
state diffusion model by extending the stochastic differential equation which governs 
the time-evolution of the wave function to the doubled Hilbert space. This procedure 
is in complete agreement with the quantum regression theorem. However, we have 
also shown that the latter fact is not sufficient to ensure that a particular simulation 
algorithm is of practical use: Although the algorithm proposed in Ref. |l3| is in 
accordance with the quantum regression theorem, it seems not to be stable in general. 
On the other hand, the scheme we proposed in this letter completely relics on the 
numerical stability of the quantum state diffusion model. 
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Figure 1. Calculation of Heisenberg operator matrix element (<j>o | <r+ (t ) | </>o ) : 
analytical solution (thick line), numerical solution using the quantum state 
diffusion unraveling of the extended quantum master equation for 10 3 realizations 
(diamonds), and the method proposed by Gisin (thin lines) for the step-sizes 
h = 0.01,0.001,0.0001. 




Figure 2. Calculation of the first order correlation function (ct+(t)<t~) s for a 
coherently driven two-level atom on resonance, (a) Analytical solution vs. the 
numerical solution (diamonds) using the quantum state diffusion model for 10 4 
realizations, (b) CPU time in seconds vs. the relative error for the simulation using 
the quantum state diffusion model (QSDE) and the quantum jump method (QJ). 
The solid lines represent the mean square deviation of the numerical solution from 
the exact solution and the dashed lines show the estimated standard deviation of 
the numerical solution. 



